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Abstract 

The time evolution of soft modes in a quantum gauge field theory is to first 
approximation classical, but the equations of motion are non-local. We show 
how they can be written in a local and Hamiltonian way in an Abelian theory, 
and that this formulation is particularly suitable for numerical simulations. 
This makes it possible to simulate numerically non-equilibrium processes such 
as the phase transition in the Abelian Higgs model and and to study, for 
instance, bubble nucleation and defect formation. This would also help to 
understand the phase transitions in more complicated gauge theories. More- 
over, we show that the existing analytical results for the time-evolution in a 
pure-gauge theory correspond to a special class of initial conditions and that 
different initial conditions can lead to qualitatively different behavior. We 
compare the results of the simulations to analytical calculations and find an 
excellent agreement. 
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I. INTRODUCTION 



The recent finding that the electroweak phase transition might be only a smooth cross- 
over |l| has changed the picture of cosmological phase transitions in a fundamental way and 
shown that theories with gauge fields can behave very differently from those with only global 
symmetries. However, the statement about the smoothness of the transition refers only to 
the equilibrium properties, and from the point of view of cosmology, it is very important to 
understand also the time evolution of the transition, in particular, whether the fields fall out 
of equilibrium and what its consequences would be. For this, a calculation scheme is needed 
for dynamical quantities that treats the gauge fields correctly and gives the same results in 
the special case of thermal equilibrium as the the dimensional reduction method that 
was used to show the smoothness of the electroweak transition. In this paper, we propose 
such a scheme. 

Instead of the electroweak theory, we will concentrate here on scalar electrodynamics, 
i.e. the Abelian Higgs model. Not only does the simpler gauge group allow more efficient 
simulations, but the theory also contains topological defects, namely Nielsen-Olesen vor- 
tices ||], and therefore it can be used to simulate defect formation in a phase transition || 
and to test the validity of the present predictions for the produced defect density || in gauge 
theories. 

In Ref. JF], the dynamics of the transition and, in particular, formation of vortices was 
simulated classically on a lattice with field equations derived from the original Lagrangian. 
A similar approach had been suggested earlier || in the context of electroweak baryogenesis. 
This was based on the observation that the quantum distribution of the long- wavelength 
modes is essentially classical. Since the time evolution of a classical field theory is given by 
the equations of motion, this approach allows non-perturbative numerical simulations. For 
baryogenesis, the interesting quantity is the hot sphaleron rate (see Ref. || and references 
therein), which is an equilibrium quantity, but similar methods can also be used to describe 
small deviations from equilibrium, such as phase transitions. After all, the phase of the 
system is only a property of the long-wavelength modes, and the distribution of the high- 
momentum modes is to a large extent independent of it. Therefore also the transition 
between the phases is described by the classical long-wavelength modes. 

This straightforward approach has the serious problem that the results depend on the 
ultraviolet cutoff |T0|JTT[] . Because of the Ray leigh- Jeans ultraviolet divergences, a classical 
continuum field theory cannot actually even be in thermal equilibrium - this was one of the 
reasons why quantum mechanics was invented in the first place! The divergences depend 
on the temperature and are absent at zero temperature, so they cannot be cancelled by 
introducing counterterms as in quantum theories. Therefore it is not sufficient to use the 
classical theory with the Lagrangian of the original quantum theory, but one has to construct 
an effective Lagrangian for the low-momentum (soft) modes by integrating out the high- 
momentum (hard) modes. The effective Lagrangian then contains precisely the necessary 
(temperature-dependent) counterterms to remove the divergences. 

In calculating static properties, i.e. expectation values of products of operators measured 
at the same time, this procedure is known as dimensional reduction ||, and leads to a three- 
dimensional effective theory, in which the temporal component of the gauge field has a Debye 
mass term m D ~ gT. To calculate non-static quantities, the full four-dimensional effective 
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Lagrangian needs to be constructed by calculating the so-called hard thermal loops, and 
in gauge theories it turns out to be non-local [p~2|, 13 . This makes computer simulations 



practically impossible, since one would have to keep in the memory all the configurations 
encountered during the time evolution. 

Fortunately, the equations of motion can be written in a local form by introducing 



auxiliary fields [0,0 • The Abelian version of the resulting system of differential equations 
is known in plasma physics as the Vlasov equation and describes collisionless electron plasma. 
Numerical solution of the equations of motion is still difficult, since the field W depends 
not only on the coordinate but also on the momenta of the hard particles it describes, 
and therefore the theory becomes essentially five- dimensional. The traditional method is 
to approximate it by a large number of point particles. This approach has been used to 
determine the hot sphaleron rate in the electroweak theory in Ref . |16[ . 

In this paper, we show how the theory can be formulated in such a way that simulations 
in terms of fields rather than particles becomes feasible. One of the problems we have to 
face is the multiplication of the number of degrees of freedom when the hard modes are 
included. We show that the number of extra degrees of freedom per lattice site must be 
chosen carefully if the simulation of a particular soft mode is to be trusted: if one wants 
to reproduce Landau damping in a mode of momentum k, naive methods require of order 
(kt) 2 degrees of freedom per lattice site for the hard modes. However, we are able to rewrite 
the effective hard thermal loop improved classical equations of motion in such a way that 
we only need of order kt extra degrees of freedom per lattice site. Our numerical methods 
reproduce known analytic results in pure Abelian gauge theory extremely well, which gives 
us confidence that we will be able to tackle the Abelian Higgs model accurately. 

The structure of the paper is the following. In Sec. || we review the hard thermal 
loop Lagrangian and derive the Hamiltonian formulation for the hard modes. In Sec. |TTI| 
we explain how to approximate the system with a finite number of degrees of freedom in 
numerical simulations. In Sec. [TV] we derive some analytical results, which we compare to 
our numerical results in Sec. [V|. Extending our technique to include the soft Higgs field 



is discussed in Sec. |\T], and Sec. |VII| contains the conclusions. The paper also has one 
appendix, in which the continuum and lattice equations of motion for the Legendre modes 
are given explicitly. 



II. KINETIC FORMULATION 

Let us consider scalar electrodynamics at a temperature T. The Lagrangian of the theory 

is 

£ = -\F, V F^ + |ZVf - m 2 |0| 2 - A|0| 4 , (1) 

where = + ieA^. We assume that e C 1 and A ~ e 2 and that m < T so that 
high-temperature approximation can be used. 

The phase structure of the theory was determined in Refs. fT7| , |TB|| . When the Higgs 
self-coupling is small, there is a first-order phase transition, but if the self-coupling is large 
enough, the transition becomes continuous. Unlike in the electroweak theory, it does not 
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become a smooth crossover. In Ref. fL9fl , it was pointed out that this can be interpreted as 
a consequence of the existence of Nielsen-Olesen vortices. 

To one-loop order, the only non- vanishing contributions to the effective Lagrangian arise 
from the two-point diagrams. The correct procedure would be to calculate the diagrams in 
the original quantum theory and to write down a classical Lagrangian that gives the same 
result when the diagrams are calculated on a lattice [[l(J. However, this is not possible in 
practice, since the form of the necessary effective lattice Lagrangian is not known and is 
presumably very complicated because the lattice has much less symmetry than the contin- 
uum. In sphaleron rate calculations, the effective Lagrangian was taken from the continuum 



theory and the parameters were fixed by matching only the static quantities [IS]. For the 
simulations presented in this paper, this problem does not arise, as will be discussed later. 
In the high-temperature approximation, the one-loop scalar self-energy is simply 

(2) 




The photon self-energy |2GJ is more complicated and if the external momentum is K = (ou, k), 
it can be written in the form 

W v = n T P£ u + IIlPl U , (3) 
where Pj? and P£ v are the spatially transverse and longitudinal projection operators 

p» = pf = o, pr^g^-^h-Pf, (4) 



and 




w , u + k s 
1 -In 



Ht = ~ (ml ~ fI L ) . (5) 



2k u-kj 2 

The Debye mass has the value m 2 D = |e 2 T 2 + 5m 2 D , where 5m 2 D is a cutoff-dependent 
counterterm. 

The self-energies can be resummed to a simple effective Lagrangian for the soft modes 



+ |£>^| 2 -,4|^| 2 -A|fl 1 , (6) 

where m 2 - = m 2 + (e 2 /4 + \/3)T 2 + Sm^, and the integration is taken over the unit sphere of 
velocities v = (l,v), v 2 = 1. Note that the form of m\ justifies using the high-temperature 
approximation even slightly below the transition. At tree level, the transition takes place 
when m\ = 0, which shows that m 2 ~ —e 2 T 2 . The minimum of the potential is therefore 
at v ~ T c , and the mass given by the Higgs mechanism is m 2 H ~ e 2 T 2 <^ T 2 . 
The equations of motion derived from the Lagrangian (^) are 

^P^U 1 - 2elm«f (7a) 
4-7T v ■ a 

D^<p = -m 2 T 4> - 2\(<p*<p)<p. (7b) 
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The derivative in the denominator in Eq. (|7a|) means that in order to simulate it numerically, 
one needs to keep the whole time evolution of the system in the memory at the same time, 
which is impossible. The form of the equation of motion (ITH) of the scalar field is much 
simpler, and actually the only effect of the hard modes is the modified mass term. The 
scalar field can therefore be simulated with standard methods and we will neglect it from 
now on, concentrating only on the gauge field. We will refer to this theory as the pure gauge 
theory although it contains the contribution from the hard scalar modes. 



The non-locality problem of Eq. ( j7a|) can be solved by introducing new fields fll4||T5| . We 
follow Ref. [jnj and add the field W(x, v), which satisfies the equation of motion 

(v-d)W(x,v) = v- E(x), (8) 
d,F^=f w (x)+f(x), (9) 



and replacing Eq. ([Ta|) with 
where 



J^W(x,v), (10) 



] w {x) = m D 



and j v {x) denotes the current due to the scalar field and any external currents. 

The system of equations (§), (H) is a special case of what is known as the Vlasov equation 
in plasma physics, where it describes collisionless electron plasma. The field W(x,v) gives 
the deviation of the density of hard particles of velocity v from their equilibrium distribution. 

In addition to the equations of motion, we also have to specify the initial conditions, and 
since we would like to have the system initially in a thermal equilibrium, they should be given 
by the equilibrium distribution, i.e. we should have a large number of initial configurations 
with a Boltzmann distribution oc exp(— /3H), where H is the Hamiltonian. However, the 
equation of motion (181) for the W field is not of canonical form, which makes it difficult to 



find the correct Hamiltonian. The suggestion of [fJTJ was to use the conserved quantity 



H= l -JSx\E 2 + B 2 + m 2 D J^W(x,v) 2 Y (11) 

In the simulations, one has to approximate the system with a finite number of degrees 
of freedom. Normally, this is done by discretizing the space to a finite lattice, but in our 
case, the field W also depends on a two-dimensional continuous parameter v. There are 
various ways to handle it in simulations. In Refs. [|l6|^], W was approximated with a 
large number of particles. The other straightforward options are approximating the velocity 
integral with a discrete set of velocities on the unit sphere, or expanding W in terms of 
spherical harmonics. However, as we will show next, at least in the Abelian case it pays to 
simplify the problem a bit first. 

Let us define in Fourier space 

f(u,k,z) = im D J ^ - / d ^ ' C W(u,k,v) (8(v- k - z) + 8(v ■ k + z 




k 2 V 2z 2 

9(u,k,z) =im D { —^W(u,k,v){5{v-k-z)+5{v-k + z)) + l ^k-A. (12) 
J A'Kv-k v ' k 2 
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Then the equations of motion in the temporal gauge Aq = are 



d 2 f(z) = z 2 V 2 f + m D z 



-V x A 



dld(z) = z 2 V • ( V0 - m D A 



8qA = -V x V x A + m D J dzz 2 ( V6» - m D A + 



1 2 U„ * /l-^ 2 



2z 2 



Vx/ + j. 



(13a) 
(13b) 

(13c) 



One advantage of this formulation is that these equations of motion are canonical and the 
corresponding Hamiltonian is 



H=\J d 3 x J dz 



E 2 + (V x A) 2 + F 2 + n 2 + z 2 {V x ff + z 2 {V9 - m D Af 



I - z 2 - - 
-m^W — - — / • V x A 



(14) 



where F = dof and II = d$9 are the canonical momenta of / and 6, respectively. We also 
need two extra conditions, namely the transverseness of / and Gauss's law 



V-/ = V-F = 0, 

V • E = -m D f dzU{z). 
Jo 



(15) 



By reformulating the degrees of freedom, we have gained two important advantages. 
First, / depends only on one internal coordinate whereas W depends on two, and second, 
since we know the Hamiltonian ( fl4]) , we know the correct distribution of the initial config- 
urations. 



III. SIMULATIONS 

We still have to approximate the ^-dependence of the hard modes with a finite number 
of degrees of freedom, and we will use Legendre polynomials for that. We define 

f n) = J o 1 dzz ] j ] ^P 2n (z)f(z), 

0V>= f 1 dzz 2 P 2n (z)9(z), (16) 
Jo 

where P n {z) is the nth Legendre polynomial. The equations of motion for different Legendre 
modes are given in Eq. ( |Al| ) in the appendix. It is also shown in the appendix that the 
approximation can be trusted if the simulation time is less than 

to « 2N max /k. (17) 

— # — * 

On the lattice, the field 9 is defined on lattice sites, while / and the gauge field A are 
defined on links between the sites in such a way that the invariance under time-independent 
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gauge transformations is preserved. Note that the field 9^ is not gauge invariant by itself, 
but the combination V6"- - 1 — ^m^A is. The canonical momenta are defined at half-way 
between the time steps so that the value of, say, / at time t + St is determined from the 
fields at time t and the momenta at time t + St/2. In this way, the time reversal invariance 
of the continuum theory is preserved. The lattice equations of motion are given in Eq. (|A4|) 
in the appendix. 

The aim of the simulations described here was mainly to test the simulation code by 
comparing it with known analytical results. Therefore we neither averaged over thermal 
initial conditions nor included the scalar field. The equations of motion for the pure gauge 
theory are linear and can therefore in principle be solved analytically. It also means that 
the modes with different momenta do not interact with each other, and to simulate a mode 
with momentum along, say, the z axis, we can use a lattice with 1 x 1 x N z sites. This 
makes the simulations very simple. Moreover, the pure gauge theory is ultraviolet finite, so 
no mass counterterms are needed. In particular, Sm 2 D and Sm 2 -, introduced in Sec. |TJ vanish. 



IV. ANALYTICAL RESULTS 

Before specializing to particular solutions, we first discuss the finite temperature prop- 
agator for the pure Abelian theory. Due to linearity, it is sufficient to consider a single 
Fourier mode at a time. Let the momentum be k, and, for simplicity, let us consider just 
the transverse part of the propagator G(t, x). Its Fourier transform has a non-analytic form 



1 



G{u,k) 



-uj 2 + k + 



m 



D 



u 2 ujk 2 



or , uj + k 



In 



UJ 



(18) 



The analytic structure of the propagator is shown in Fig. |l|. The propagator has two os- 
cillatory poles on the real axis at uj 1 



k 2 + \m 2 j 



t D , but also a branch cut from uj = — k to 
uj = k. The original equation of motion ([7a] ) implies that the branch cut must be taken 
along the real axis. There are no other singularities on this physical Riemann sheet, but 
if one continues the propagator analytically from the upper half-plane through the branch 
cut, one finds a pole at 



4k 3 



U = -111 



Tim 



2 ' 

D 



(19) 



Likewise, continuing analytically from the lower half-plane through the cut reveals a pole at 
uj = ij L . 

The solution to the inhomogeneous equation with an external source j(uj) is, assuming 
that the current is transverse and that the fields vanish at t = — oo, 



A(t, k) 



oc+ie du) 
-oo+ie 2,7T 



e-^G(u)j(uj). 



(20) 



If t > 0, the integral can be transformed into a contour integral by closing it on the lower 
half-plane. This integral has three pieces: two contributions from the poles and one from 
the integral around the branch cut. The result of the integration around the cut is sensitive 
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to the functional form of a reflection of the non- locality of Eq. (£a|): the behavior of 

the system depends on its whole history. 

For example, Boyanovsky et al. |23j discussed a situation in which an inhomogeneous 
initial configuration is set up for the soft fields A and E, and calculated its relaxation to the 
equilibrium at asymptotically long times. Still keeping the assumption of transverseness, 
the solution was given as a function of A and E at time t = 0, 



A(t, k) « A(0,k) 
where 



4 cos kt 

Z[T\ cos u p t = -— 

m D t 2 



E(0,k) 



grj,. sin u p t 4 sin kt 
lj p m 2 D kt 2 



(21) 



z ^=-(^y • <22) 

Eq. (pi| ) can be seen to correspond to j(uS) = — iuA(0) — E(0). As they emphasize, the 
dominant contribution does not come from the smallest frequencies but those near the end 
points of the branch cut. 

With a different choice of j, qualitatively different solutions can be found. Suppose that 
j(t) is increased very slowly from zero to a finite value and then suddenly switched off, 
i.e. j(t) = j e lot Q(—t), where 70 is eventually taken to zero and 0(t) is the step function. 
The Fourier transform j(uS) = jo/(lo + iw) is peaked around the origin, and therefore the 
integral around the branch cut is dominated by the integrand near uj = 0. Hence, expanding 
the propagator around uj = ie and u = —ie on the upper and lower half-planes, respectively, 
we may write 

Mt, k)\ t « Jo f p. ( e -~-L — L_ + **_L. — LJ) (23) 

cut J-k in \ 70 + iuj 7l — ioj 70 — ioj 7l — ioj J k 

The limits of the integral may be taken to infinity, and the contour closed in either the 
upper or the lower half-plane depending on the sign of the exponent. The second term has 
no poles in the upper half-plane and hence vanishes, while the first term gets a contribution 
from the pole at uj = — iji. Thus we find an exponentially damping solution 

X ^L M -^'- < 24 > 

This exponential decay is known as Landau damping, and it is important to notice that it 



did not really come from integrating around a pole of the propagator dig) . If it had, then 
the equations of motion would either have to break time reversal invariance, which they 
do not do, or there would have to be a corresponding exponentially growing solution, thus 
making the system unstable. 

For completeness, we also calculate the contribution coming from the poles at uj = ±uj p . 
The full result is then 

At, k) « -Jo (z[Tf-^ + le-^) . (25) 



S 



V. NUMERICAL RESULTS 



In order to numerically reproduce the results fl2~ID , ([2~5"D, we have to specify the correct 
initial conditions in terms of / and 9. The result ([H]) is given simply by f(0) = 6(0) = 0. 
For the soft modes, we chose A(0) = 0, E(0) = 10?/ sin A; • x. To see the power-law damping 
most clearly, mjj must be relatively small, and we chose = 2tt, k = 2ttx in our units. 
The lattice size was 20 x 1 x 1, lattice spacing a = 0.05 and the time step 5t = 0.01. We 
used A max = 200 so that according to Eq. (|17]) the results should be reliable when t < 60. 
Eq. ( pip becomes 

A{t, k) « (-1.236 sin 7.451x + 0.1613 ^^ ) y. (26) 

We added to the numerical result a constant-amplitude part a sin (3t and determined the 
parameters a and /3 from the condition that the remaining part agrees with the decaying 
part of Eq. (HI). The result with a = 1.24008 and (3 = 7.43024 is shown in Fig. |. The 
discrepancy between the numerical and analytical results is less than one percent. 

The initial conditions that correspond to Eq. (|25f) are those in which the first and second 
time derivatives of the fields vanish. If we choose jo such that A(0,x) = ysink ■ x, then 
Eqs. fll|) and fllBD imply that 



(0, x) = -zy-z cos k ■ x, 



3k 
m D 



(0, x) = ^-z cos k ■ x, 

15k 



\n>l) 



(0,x) = 0, ^^(0,f)=0. (27) 



The exponential damping is seen most clearly if 3> k, and we chose two different values 
m£> = 107r and = 207r while we took again k = 2ttx. Again, the lattice size was 20 x 1 x 1, 
the lattice spacing a = 0.05, the time step 5t = 0.01, and A max = 200. The results of the 
simulations are shown in Fig. |3|. The predicted decay rates are jl = 0.32 and jl = 0.08 and 
the amplitudes of the oscillation A osc w 0.103 and A osc m 0.029. As can be seen from Fig. |^, 
these agree very well with the numerical results. 

We also carried out simulations with with different values of A max to find out at what 
time the approximation breaks down and to test the estimate (|TTD. We used mp = 20tc 
and the values for the other parameters were as before. Since all the Legendre modes with 
n > 1 are initially zero, we expect the first errors to occur at time t ~ 2t w AN mSuX /k. The 
result is shown in Fig. f| for various values of iV max and confirms our estimate. The number 
of extra scalar degrees of freedom needed for simulating with any particular value of iV max 
is 8A max . 

In order to compare the efficiency of the Legendre polynomial formulation with other 
approaches, we carried out the same simulations with a more straightforward method. We 
chose a large number of points on the unit sphere to represent different velocities and used 
them to simulate the pair of equations (|8|), @. More precisely, the different velocities were 

(N v + i, n 1n n z ) , , 

' ' " -N v < m < N v , (28) 



(N v + \f + n 2 y + n 
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and those obtained from that with reflections or rotations of ir/2. The parameters were 
the same as in Fig. ^. The values of N v used ranged from 2 to 8, and the corresponding 
number of extra degrees of freedom is 6(2N V + l) 2 . The results in Fig. [| show clearly that 
the number of different velocities becomes quickly prohibitive when the simulation time is 
increased. If we assume that, as with the Legendre polynomials, the particularly smooth 
initial conditions result in a factor of two in the reliable simulation time, we can estimate 
that generally a simulation will be reliable if t ^ irN v /k. 

This estimate can also be reached analytically. Since the time evolution is a superposition 
of oscillations, the reliable simulation time is t « tt/ujq, where Uq is the smallest frequency, 
i.e. the pole of the propagator that is nearest to the origin. The transverse self-energy 
Ht(co>) diverges if u — k ■ v for any velocity v. The smallest pole of the propagator is 
at a frequency ujq, which is smaller than min{fc • v}, but of the same order of magnitude. 
Generally min{&; • v} « k/N v , leading to the previous estimate. 

Since expansion of W in Eq. (|8]) in terms of spherical harmonics Y™ is essentially equiv- 
alent to the expansion of / and 9 in terms of Legendre polynomials, we can also estimate 
the effectiveness of that approach. In order to reach the same simulation time, the highest n 
used should be 2iV max . The number of extra degrees of freedom would then be (2iV max + l) 2 , 
which again increases much faster than in the Legendre polynomial approach. Our estimates 
of the reliable simulation times in the different approaches have been summarized in table |. 
They show that the Legendre polynomial formulation is the most economical one. [] 



VI. HIGGS MODEL 

The pure gauge theory that has been discussed so far is an ideal way to test the formalism, 
since it is linear and one can in principle solve it exactly. In the previous section, we have 
shown that the numerical results agree with the analytical calculations and we were even 
able to estimate the number iV max of Legendre modes needed to ensure that the simulation 
is reliable. 

It is very easy to transform the system to a very non-trivial one simply by adding a 
Higgs field. While this makes it impossible to solve the equations of motion analytically, 
from the point of view of numerical simulations it means only adding one more field. Still, 
the resulting theory has a complicated structure with a phase transition and topological 
defects. 

As was shown in Sec. Q, the weak-coupling condition e C 1 and A ~ e 2 implies that the 
hard thermal loop approximation is also valid near the transition in the broken phase. The 
reason is that the distribution of the hard modes is independent of the phase. Furthermore, 



1 The analogous number of degrees of freedom in the particle method used by Moore et al. [16] 
is 6(n). The value of (n) they were using varied between 17 and 120, but they did not carry 
out this kind of systematic analysis of the corresponding reliable simulation time. However, the 
intrinsic randomness of the method seems to imply that it does not reproduce the correct behavior 
as accurately at short times as the methods discussed here. 
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when the phase transition takes place in a finite time, only the soft modes fall out of equi- 
librium, and therefore the classical description remains valid, even in such non-equilibrium 
processes. 

The above suggests that the classical formalism discussed in this paper can be used to 
simulate non-equilibrium dynamics of the phase transition in hot scalar electrodynamics. In 
the cosmological setting the transition takes place when the universe cools as it expands, 
but in a simulation some other way to change the temperature is needed. A straightforward 
way of doing this is preparing an initial ensemble with a higher temperature for the soft 
modes than for the hard modes. A large number of initial configurations would then be 
taken from this ensemble and evolved in time. When the soft and hard modes interact, the 
temperature of the soft modes decreases and they undergo a phase transition. 

While this kind of an instantaneous quench in the temperature is not very realistic, it 
would still be a first approximation to the true cosmological phase transition. These simula- 
tions would answer many interesting questions about the dynamics of the phase transition, 
for instance the density of topological defects created. 

Certainly, the hard thermal loop approximation we have used is only valid at relatively 
short times. We have neglected the collisions of the hard particles and at long times t ^ 
l/e 4 T they would have an essential contribution to the time evolution If the coupling 
is weak, as we assume, this is not a serious problem. The other problem is to remove the 
effect of the ultraviolet lattice modes as discussed in Sec. |I|. The fact that the simulations 
of [|7| showed no signs of a Rayleigh- Jeans catastrophe developing over the course of the 
simulation seems to indicate that the effect can safely be neglected. 

Phase transitions and other non-equilibrium processes can be simulated also in non- 
Abelian theories using the Vlasov equation (|[) to describe the hard modes. It would be 
important to know, whether a formulation analogous to Eq. ([Tj^) with only one internal 
coordinate is also possible in non-Abelian theories, since that would reduce the need of 
computational power drastically. Since the derivative in Eq. (|8]) is replaced with a covariant 
derivative, operating in the momentum space becomes more complicated and the same 
derivation cannot be used. 



VII. CONCLUSIONS 

In this paper, we have studied hot Abelian gauge field theory in the hard thermal loop 
approximation. Starting from the local kinetic formulation (|8p, ([|), we were able to refor- 
mulate the degrees of freedom in a more economical way. The equations of motion in the 



new formulation are canonical and we found the explicit form of the Hamiltonian flTl]). 

The pure gauge theory discussed in this paper is linear and can therefore be solved 
analytically, at least in principle. We pointed out that because of non-locality it is not 
sufficient to specify the initial conditions for the soft modes only. In a sense, one will have 
to specify the whole history of the system in order to calculate its behavior in the future. 



Taking this effect into account modifies the existing results [23[ and leads to exponential 
damping (^) with the Landau damping rate. 

In order to simulate the system, we approximated the functions /, 9 with a finite number 
of Legendre polynomials. The simulations reproduced the analytical results, and the number 
of degrees of freedom needed to describe the hard modes was much smaller than in other 
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possible approaches. This shows that the formulation presented in this paper is very well 
suited for numerical simulations. 

While the pure gauge theory is trivial in the sense that it can be solved analytically, 
including the soft Higgs field makes analytical calculations essentially impossible. On the 
other hand, it is very simple to add it to numerical simulations. Leaving possible ultraviolet 
problems aside, that would give a way to study non-perturbatively the non-equilibrium 
dynamics of a phase transition in a gauge field theory. 
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APPENDIX A: EQUATIONS OF MOTION FOR LEGENDRE MODES 

The equations of motion for the Legendre modes defined in Eq. fll6|) are 

n-l) 



A t( 1 c 1 c 4 c 

+m D V X A —0 nfi + -rrzOn,! - 7TiT^n,2 

V15 105 315 

3*0(») = c+v 2 # (n+1) + c°v 2 # (n) + c-v 2 e (ft -^ 

^ 7( l c 4 c 8 c 

-m D V ■ A -d n „ + ~<Vi + ^77^,2 
V5 35 Sib ' 



9 2 I=-VxVxI - ^m 2 D A + m D (w {0) + V x f 0) ) , (Al) 



where 



c+ _ (2n+l)(2n + 2) ^ _ 1 / (2n + l) 2 | 4n 2 \ Q _ 2n(2n - 1) 



(4n+l)(4n + 3) n An + 1 \ An + 3 An - 1 J ' n (4n + l)(4n - 1) 

(A2) 

Approximating an infinitely many degrees of freedom with a finite number gives neces- 
sarily rise to errors. Suppose we only take N mSu ^ lowest Legendre modes into account. If 
n>l, then the equation of motion for the Fourier mode with momentum k of 9^ is simply 

d$e (n) = ~k 2 (# (n+1) + 20< n > + # (n - 1} ) . (A3) 

By writing 6^ = {—l) n 6^ n \ we notice that this is precisely a discretized version of a wave 
equation, with waves propagating at speed c — k/2. The same applies to / as well. Since in 
the approximation, errors only occur at n — N max , we can estimate that the approximation 
works as long as the error has not had time to propagate to n — 0, which is the only mode 
that is coupled to the observable soft modes. Thus, for a mode with momentum k, the 
approximation is reliable as long as the simulated time is smaller than 
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to « 2N nmx /k. (|T3) 

Therefore, if we want to measure correlators with time separation t max , we need to have 
N max ~ t max /2a, where a is the lattice spacing. 

For completeness, we give here the lattice equations of motion 

St ^ / St _ 

^i(t+2-,a;) = ^i(t-2-,x) 

+5t {^AT^-(t, f) + im^,(t, £) - — [A t +0(°> + e, jfc AT/i 0) (t, x) 
l„ a o a 

F i (t + —,x) = F> '{t- — ,x) 

+St! [ C+V*ft +1) (t,x) + C^ft\t,x) + C-V 2 ft 1] (t,x) 
m D ( 1 1 4 



m,D ( 1 r 1 r 4 \ . , . , . 

+— (^n,o + — ^ - — 5„, 2 j A fc (t,s) 



nW(t + |,f) = nW(t-|,f) 



+^jc+v 2 # (n+1) (M) + c°v 2 # (n) (M) + ^v 2 ^ (n - 1) (t, 

m D (\ 4 r 8 . \ . _ . , „ 1 

:0n,0 + ~On,l + T^°n,2 A { Aj(t, X) S, 



x 



a \5 ' 35 ' 315 
+ <5t, x) = Ai(t, x) - 5tEi(t 

ft\t + St, x) = ft\t, x) + 6tFt\t + |, x), 

9 (n) (t + St,x) =6 {n \t,x) + 5m( n \t+j,x), (A4) 

where we have used shorthand notations 

A±0(x) = ±(0(x±*)-0(x)), 
A i:j (x) = AfAj(x) - A+Mx), 

V 2 0(x) = ^ E f + *) - 2( f>(^ * ) + x - z)) . (A5) 
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FIGURES 




Re uj 



FIG. 1. Analytic structure of the propagator (18). The filled circles are oscillatory poles, the 
thick line is the branch cut, and the arrows show how imaginary poles at uj = ±ijL, depicted by 
open circles, can be found by analytically continuing the propagator through the branch cut. 
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FIG. 2. The comparison of Eq. (gray line) to the numerical result (black line) with 

k = mo = 27T. The constant-amplitude oscillation has been subtracted. 
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FIG. 3. Results of the Landau damping simulations with k = 2tt and = Wtt (left) and 
mo = 207T (right). The white line is the predicted Landau damping rate jl and the dashed lines 
show the envelope of the analytical result (25). 
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FIG. 4. Landau damping simulations with iV max = 10 . . . 40 (k = 2tt, mj) = 207r). The vertical 
dashed lines show the analytical prediction t ~ 4A r max /fe for the time when the approximation 
should break down. The black curve is the result with -/V max = 200 and the gray curves correspond 
to iV max = 10, 20, 30, 40 from left to right. 
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FIG. 5. Landau damping simulations with k = 2ir, mo = 20tt in a formulation in which a large 
number of points on the unit sphere of velocities is used to approximate the integral. The values 
of iV„ were (from left to right) 2, 4, 6 and 8 and the corresponding numbers of degrees of freedom 
are shown in the plot. 
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TABLES 



TABLE I. The approximate number of degrees of freedom needed in various approaches to 
simulate a mode with momentum k reliably for time t S> 1/k. 



Approach 


Degrees of freedom 


Legendre polynomials 


4kt 


Spherical harmonics 


(ktf 


Discrete velocities 
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